The problem on 2-dimensional grid is more interesting that 1-dimensional problem. In this case direct substitution (as in 1-dimensional case) doesn’t work.
If all coefficients σ = (σj)j=0,…,n, in matrix A2 are equal, then there are at least two ways to solve the system of equations using nearly О(N) floating point operations, where N = (n + 1)2 – the size of the matrix in 2-dimensional case.
The first way is to use fast Fourier transform. Indeed, when σ is constant, matrix А2 is a Toeplitz matrix (discrete variant of convolution). After applying Fourier transform, convolution becomes an operator of multiplication on a function. This operator is easy to invert. After that it remains to apply inverse Fourier transform.
As soon as fast Fourier transform is a fast operation, matrix А2 may be inverted using O(N log(N)) operation. In other words, we have nearly О(N) complexity (i.e. nearly the same degree of complexity as in 1-dimensional case).
Another way to solve system of equations with matrix А2 is to do nearly the same as we did in 1-dimensional case – use direct substitution. The only difference is that we should do corresponding arithmetic operations (add, multiply, divide) not with numbers but with blocks (these blocks in 2-dimensional problem are very simple matrices). This approach is called cyclic reduction.
One more thing that helps here, is that all blocks (when σ is constant) on the main diagonal of the matrix А2 are the same. Thus, we need to invert corresponding matrix of a block only once.
Accordingly, for constant σ we can solve the system of equations with matrix А2 for О(N) operations (i.e. complexity of the problem is O(n2)).